Import data

# Import data
df <- readr::read_rds('./data/Juliane_raw_FEST.rds')

Quick look

# Quick look
dim(df)
## [1] 2735   11
head(df)
## # A tibble: 6 x 11
## # Groups:   id [1]
##      id block_number trial_no   trial_type perc_trial_type      CS    US
##   <dbl>        <chr>    <int>        <chr>           <chr>   <chr> <chr>
## 1     8     baseline        1      CS_only         CS_only  CSplus  <NA>
## 2     8     baseline        2 CSminus_only    CSminus_only CSminus  <NA>
## 3     8     baseline        3      CS_only         CS_only  CSplus  <NA>
## 4     8     baseline        4      CS_only         CS_only  CSplus  <NA>
## 5     8     baseline        5 CSminus_only    CSminus_only CSminus  <NA>
## 6     8     baseline        6 CSminus_only    CSminus_only CSminus  <NA>
## # ... with 4 more variables: rating <int>, rating_time <int>,
## #   button_correct <chr>, button_time <int>
tail(df)
## # A tibble: 6 x 11
## # Groups:   id [1]
##      id block_number trial_no    trial_type perc_trial_type      CS    US
##   <dbl>        <chr>    <int>         <chr>           <chr>   <chr> <chr>
## 1    30        test3      137 CSminus_USlow   CSminus_USlow CSminus   low
## 2    30        test3      138  CSminus_only    CSminus_only CSminus  <NA>
## 3    30        test3      139 CSplus_UShigh   CSplus_UShigh  CSplus  high
## 4    30        test3      140  CSminus_only    CSminus_only CSminus  <NA>
## 5    30        test3      141 CSminus_USlow   CSminus_USlow CSminus   low
## 6    30        test3      142 CSplus_UShigh   CSplus_UShigh  CSplus  high
## # ... with 4 more variables: rating <int>, rating_time <int>,
## #   button_correct <chr>, button_time <int>
glimpse(df)
## Observations: 2,735
## Variables: 11
## $ id              <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8...
## $ block_number    <chr> "baseline", "baseline", "baseline", "baseline"...
## $ trial_no        <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,...
## $ trial_type      <chr> "CS_only", "CSminus_only", "CS_only", "CS_only...
## $ perc_trial_type <chr> "CS_only", "CSminus_only", "CS_only", "CS_only...
## $ CS              <chr> "CSplus", "CSminus", "CSplus", "CSplus", "CSmi...
## $ US              <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "low",...
## $ rating          <int> 0, -7, -7, -24, -32, -39, -17, -33, -40, -33, ...
## $ rating_time     <int> 4824, 4342, 5833, 4173, 2689, 3275, 3460, 3201...
## $ button_correct  <chr> "correct", "correct", "correct", "correct", "c...
## $ button_time     <int> 8407, 5485, 7576, 5572, 3792, 4610, 4875, 4152...

Plot ratings by trial type

Ratings on the FEST are plotted for each participant and then for the whole sample. First we plot by trial type delivered, then by trial type perceived.

# Make plots
ggplot(data = df) +
  aes(x = trial_no,
      y = rating,
      colour = trial_type) +
  geom_point() +
  scale_colour_brewer(palette = "Set1") +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 10.5, linetype = 3) +
  geom_vline(xintercept = 34.5, linetype = 3) +
  facet_wrap(~ id) +
  labs(title = "Ratings over course of procedure",
       subtitle = "Faceted by ID")

df_test <- df %>% filter(trial_type == "CSplus_UStest" | trial_type == "CSminus_UStest")

ggplot(data = df_test) +
  aes(x = trial_no,
      y = rating,
      colour = trial_type) +
  geom_point() +
  scale_colour_brewer(palette = "Set1") +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 70.5, linetype = 3) +
  geom_vline(xintercept = 34.5, linetype = 3) +
  geom_vline(xintercept = 105.5, linetype = 3) +
  facet_wrap(~ id) +
  labs(title = "Ratings over course of test phases, UStest trials only",
       subtitle = "Faceted by ID")

ggplot(data = df_test) +
  aes(x = as.factor(trial_type),
      y = rating,
      colour = trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  facet_wrap(~ id) +
  labs(title = "Ratings of test stimuli according to actual trial type",
       subtitle = "Faceted by ID")

ggplot(data = df_test) +
  aes(x = as.factor(perc_trial_type),
      y = rating,
      colour = perc_trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  facet_wrap(~ id) +
  labs(title = "Ratings of test stimuli according to perceived trial type",
       subtitle = "Faceted by ID")

df_test %>% filter(!is.na(trial_type)) %>%
ggplot(data = .) +
  aes(x = as.factor(trial_type),
      y = rating,
      colour = trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  labs(title = "Ratings of test stimuli",
       subtitle = "According to actual trial type")

df_test %>% filter(!is.na(perc_trial_type)) %>%
  ggplot(data = .) +
  aes(x = as.factor(perc_trial_type),
      y = rating,
      colour = perc_trial_type) +
  geom_jitter(height = 0) +
  geom_boxplot() +
  scale_colour_brewer(palette = "Set1") +
  labs(title = "Ratings of test stimuli",
       subtitle = "According to perceived trial type")

Prediction interval plots

Here, we plot the prediction interval for the ‘control’ stimulua type (e.g. CSminus/UStest) and depict the trials of the comparator trial type (e.g. CSplus/UStest) that fall within and outside the predictio interval. We do this for the test trials, as delivered and as perceived, and then for the acquisition trials.

# Prediction intervals for test trials as DELIVERED
df_prediction <- df_test %>%
  # Filter for CSminus_UStest
  filter(trial_type == 'CSminus_UStest' & id != 23) %>%
  # Calculate 95% prediction interval 
  group_by(id) %>%
  dplyr::summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
            upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
  # Truncate PIs at +50 and -50
  mutate(lower_pi = ifelse(lower_pi < -50,
                           yes = -50,
                           no = lower_pi),
         upper_pi = ifelse(upper_pi > 50,
                           yes = 50,
                           no = upper_pi)) %>%
  # Join back with df_test
  left_join(df_test) %>%
  # Is CSplus_UStest within the prediction interval
  mutate(CC_effect = ifelse(rating > upper_pi,
                               yes = 'yes',
                               no = 'no')) %>%
  # Filter for CSplus_UStest
  filter(trial_type == 'CSplus_UStest') 

ggplot(data = df_prediction) +
  aes(x = as.factor(id),
      y = rating) +
  geom_jitter(aes(colour = CC_effect),
             height = 0) +
  geom_errorbar(aes(ymin = lower_pi,
                    ymax = upper_pi)) +
  ylim(-50, 50) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_colour_brewer(palette = "Set1") +
  coord_flip() +
  labs(title = 'Prediction interval for CSminus_UStest trial ratings',
       subtitle = 'Points are individual CSplus_UStest trial ratings',
       y = 'FEST rating',
       x = 'Participant ID')

# Prediction intervals for test trials as PERCEIVED

id_19 <- df_test %>% filter(id == 19)

plot(id_19$rating, id_19$trial_no)

df_prediction_perc <- df_test %>%
  filter(!is.na(perc_trial_type)) %>%
  # Filter for CSminus_UStest
  filter(perc_trial_type == 'CSminus_UStest' & id != 23 & id != 19) %>%
  # Calculate 95% prediction interval 
  group_by(id) %>%
  summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
            upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
  # Truncate PIs at +50 and -50
  mutate(lower_pi = ifelse(lower_pi < -50,
                           yes = -50,
                           no = lower_pi),
         upper_pi = ifelse(upper_pi > 50,
                           yes = 50,
                           no = upper_pi)) %>%
  # Join back with df_test
  left_join(df_test) %>%
  # Is CSplus_UStest within the prediction interval
  mutate(CC_effect = ifelse(rating > upper_pi,
                            yes = 'yes',
                            no = 'no')) %>%
  # Filter for CSplus_UStest
  filter(perc_trial_type == 'CSplus_UStest') 

# Save file df_prediction_perc

write_rds(x = df_prediction_perc,
          path = './data/Juliane_raw_as_perc.rds')

ggplot(data = df_prediction_perc) +
  aes(x = as.factor(id),
      y = rating) +
  geom_jitter(aes(colour = CC_effect),
              height = 0) +
  geom_errorbar(aes(ymin = lower_pi,
                    ymax = upper_pi)) +
  ylim(-50, 50) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_colour_brewer(palette = "Set1") +
  coord_flip() +
  labs(title = 'Prediction interval for ratings of trials perceived as CSminus_UStest ',
       subtitle = 'Points are individual CSplus_UStest (percieved) trial ratings',
       y = 'FEST rating',
       x = 'Participant ID')

# Prediction intervals for all reinforced (CS+/USH or CS-/USL) trials

df_prediction_reinf <- df %>%
  # Filter for CSminus_USlow
  filter(trial_type == 'CSminus_USlow' & id != 23) %>%
  # Calculate 95% prediction interval 
  group_by(id) %>%
  dplyr::summarise(lower_pi = mean(rating, na.rm = TRUE) - 1.96 * (sd(rating, na.rm = TRUE)),
            upper_pi = mean(rating, na.rm = TRUE) + 1.96 * (sd(rating, na.rm = TRUE))) %>%
  # Truncate PIs at +50 and -50
  mutate(lower_pi = ifelse(lower_pi < -50,
                           yes = -50,
                           no = lower_pi),
         upper_pi = ifelse(upper_pi > 50,
                           yes = 50,
                           no = upper_pi)) %>%
  # Join back with df_test
  left_join(df) %>%
  # Is CSplus_UShigh within the prediction interval
  mutate(diff_intensities = ifelse(rating > upper_pi,
                            yes = 'yes',
                            no = 'no')) %>%
  # Filter for CSplus_UShigh
  filter(trial_type == 'CSplus_UShigh') 

ggplot(data = df_prediction_reinf) +
  aes(x = as.factor(id),
      y = rating) +
  geom_jitter(aes(colour = diff_intensities),
              height = 0) +
  geom_errorbar(aes(ymin = lower_pi,
                    ymax = upper_pi)) +
  ylim(-50, 50) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_colour_brewer(palette = "Set1") +
  coord_flip() +
  labs(title = 'Prediction interval for CSminus_USlow trial ratings',
       subtitle = 'Points are individual CSplus_UShigh trial ratings',
       y = 'FEST rating',
       x = 'Participant ID')

Conclusion about calibration

Calibration was poor. However, participants 11, 15 and 26 have a prediction interval for their CS-/USlow trials that does not cross zero, and they also show no CC effect according to the previous plot.

Session information

sessionInfo()
## R version 3.4.0 (2017-04-21)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2       RColorBrewer_1.1-2 magrittr_1.5      
##  [4] dplyr_0.7.4        purrr_0.2.4        tidyr_0.7.2       
##  [7] tibble_1.3.4       ggplot2_2.2.1      tidyverse_1.1.1   
## [10] readr_1.1.1       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.13     cellranger_1.1.0 compiler_3.4.0   plyr_1.8.4      
##  [5] bindr_0.1        forcats_0.2.0    tools_3.4.0      digest_0.6.12   
##  [9] lubridate_1.6.0  jsonlite_1.5     evaluate_0.10.1  nlme_3.1-131    
## [13] gtable_0.2.0     lattice_0.20-35  pkgconfig_2.0.1  rlang_0.1.2     
## [17] psych_1.7.8      yaml_2.1.14      parallel_3.4.0   haven_1.1.0     
## [21] xml2_1.1.1       httr_1.3.1       stringr_1.2.0    knitr_1.17      
## [25] hms_0.3          rprojroot_1.2    grid_3.4.0       glue_1.1.1      
## [29] R6_2.2.2         readxl_1.0.0     foreign_0.8-69   rmarkdown_1.6   
## [33] modelr_0.1.1     reshape2_1.4.2   backports_1.1.1  scales_0.5.0    
## [37] htmltools_0.3.6  rvest_0.3.2      assertthat_0.2.0 mnormt_1.5-5    
## [41] colorspace_1.3-2 labeling_0.3     stringi_1.1.5    lazyeval_0.2.0  
## [45] munsell_0.4.3    broom_0.4.2